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Abstract 



A new theory for the dynamics of the magnetic particles and their mag- 
netic moments in ferrofluids is developed. Based on a generalized Lagrangian 
formulation for the equations of motion of the colloidal particle, we intro- 
duce its interaction with the solvent fluid via dissipative and random noise 
torques, as well as the interactions between the particle and its magnetic mo- 
ment, treated as an independent physical entity and characterized by three 
generalized coordinates, its two polar angles and its modulus. It has been 
recognized recently that inertial effects, as well as the particle's rotational 
Brownian motion, may play important roles on the dynamic susceptibility of 
a class of magnetic fluids. No satisfactory theory existed, up to now, that 
takes this effects into account. The theory presented here is a first-principles 
3-dimensional approach, in contrast to some phenomenological 2-dimensional 
approaches that can be found in the recent literature. It is appropriate for 
superparamagnetic, non-superparamagnetic and mixed magnetic fluids. As a 
simple application, the blocked limit (magnetic moment fixed in the particle) 
is treated numerically. The rotational trajectory of the particles in presence of 
a magnetic field, as well as the response functions and dynamic susceptibility 
matrices are explicitly calculated for some values of the parameters 
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I. INTRODUCTION 



Considerable interest has been shown, in recent years, on the dynamics of the magneti- 
zation of ferrofluids in presence of applied magnetic fields and on the corresponding complex 
magnetic susceptibility. Just to give a few examples of recently published work on the field 
we mention the theoretical works by Raikher and Rusakov [|IJ, Coffey and Kalmykov ||, 
Shliomis and Stepanov ||], the experimental works by Morais et al. 0], Fannin et al. ||, 
Vincent et al. M and an experimental-theoretical paper by Fannin et al. J7| . Certainly, this 
increased interest in a better understanding of the behavior of these materials is related to 
their renewed technological importance, with various new applications ||. 

The usual theoretical approach to calculate the dynamic susceptibility is based on 
Gilbert's || or Landau Lifshitz' |]TU[ equation (which are equivalent) for the dynamics of 
the magnetic moment, with the addition of noise, following the pioneering work of Brown 
|TIJ . Several authors use these equations of motion to calculate relaxation times and the 



susceptibility is then borrowed from Debye's theory fT2|| . 

Two distinct rotational relaxation mechanisms may coexist in ferrofluids: the Neel re- 
laxation, by which the magnetic moment moves with respect to the mechanical particle, 
and the Brownian, or Debye relaxation, corresponding to the particle's rotation inside the 
fluid. In most experimental situations one of these mechanisms is dominant, and this may 
be the reason why there is not, up to now, that we know, a satisfactory theory, sufficiently 
general to be applied for all situations, from the pure Neel to the pure Brownian relaxation, 
passing by all possible combinations of those mechanisms. In this respect the model of "two 
spheres", by Fannin and Coffey fl3| should be mentioned as a first effort. 

The purpose of the present paper is to present such a general theory. The main limitation 
of our approach is that we deal only with axially symmetric particles, with easy axis of 
magnetization parallel to the symmetry axis. However, the magnetic moment is allowed to 
rotate inside the particle, as well as to have an oscillating modulus, and the particle is allowed 
to rotate with respect to the solvent, which is immobile with respect to the laboratory. 
The suspension is considered sufficiently dilute for the particle-particle interaction to be 
negligible, so that we deal only with single particle dynamics. However, in a mean field 
approximation, our approach can serve as a starting point for the inter-particle interactions 
to be considered in future works. 

In section II we write the equations of the rotational motion of an axially symmetric 
particle inside a fluid (Langevin-type equations), based on the generalized Euler-Lagrange 
equations. In section III we obtain, from the equations of section II, in a convenient limit, 
the equations of motion for the magnetic moment /x, which reduce, in the case of constant 
modulus of n, to the Gilbert's equation. In section IV we arrive at the set of six coupled 
equations, for the six degrees of freedom, the three Euler angles of the particle's rotations, 
the two polar angles of /x and its modulus. Some less general situations are also considered 
in this section, as particular cases. In section V the "blocked limit", i.e., when the magnetic 
moment is fixed with respect to the particle, is treated as an explicit example of application. 
In section VI we introduce a simple version of linear response theory, applicable for the cases 
where the noise can be considered only for its effect as a thermal bath. In section VII we 
apply this linear response approach to calculate the dynamic susceptibility of the ferrofluid 
in the blocked limit and in section VIII some numerical results are presented and discussed. 
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We do not explore, in the present paper, the set of six equations of section IV, Eqs. (|llD, 
in its great generality, because this would make the paper too long. Work with this purpose 
is being carried out by the authors, to be published in future papers. 



II. ROTATIONAL DYNAMICS OF A PARTICLE IN A FLUID 

Consider a particle of axially symmetric shape in suspension in a fluid. The principal 
moments of inertia will be denoted by I\ = I2 and I3. Disregarding translational degrees 
of freedom, its Lagrangian may be written in terms of the Euler angles 9, <p and ip (in the 



notation of Goldstein taken as generalized coordinates, as 



L = |(# 2 + 2 sin 2 9) + + cos 9) 2 - V{9, 0) (1) 

where V(9, 0) is some orientation dependent potential. It cannot depend on ip because of 
the axial symmetry of the particle. 

The interaction forces (torques) between the particle and the fluid are of the dissipative 
and noise types. Therefore, they are not included in the Lagrangian, but instead, we have to 
use the "generalized Euler-Lagrange equations" , with the corresponding torques, represented 
by Qi, at the right hand side: 

d dL dL 

dt d<ji dqi 

where qi = 9,<p or ip. 

We write the non-conservative torques Qi as sums of dissipative and noise terms, in the 
form 

Qi = ~+T i (t) i (3) 
oqi 



where T is the following Rayleigh dissipation function [14 



T = ^\((9 2 + 2 sin 2 9) + l - A'(V» + cos 9)\ (4) 

and Ti(t) are the noise torques. The dissipation constants A and A' may be different because 
A' is associated with the particle rotation around the symmetry axis, while A is associated 
with the rotations perpendicular to it. Substituting Eqs. (]l|), @ and (P into Eq. (||) we 
obtain the following system of equations for the particle's rotation: 

h (9 - 2 sin 9 cos 9) + J 3 (ip + cos 9) sin 9 + A 9 + V e = T e , (5a) 
Ji^sin 2 ^ + 2 9 sin 9 cos 9) + J 3 cos6) — (ip + 0cos0) + 

-h(ip + cos 9)9 sin ^ + A sin 2 6 + V4 = , (5b) 

hj-(ij> + 0cos6») + A' (ip + 0cos6») = r^. (5c) 
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where Vq = dV/39 and = dV/dcj). The expression (tp+cf) cos 9) was left unbroken wherever 
it appears in the above equations because it represents the component of the angular velocity 
vector u) along the symmetry axis and we make use of this fact in the interpretation of the 
dissipative torques in terms of the components of u>, as follows. 

Let us define the following four unit vectors: z, along the laboratory z-axis, c, along the 
particle's symmetry axis, a, perpendicular to the plane containing c and z (cz-plane) and 
b, perpendicular to the ca-plane, namely, 

* = (0, 0, 1) , (6a) 
c— (sin 9 cos 0, sin 9 sin 0, cos 9) , (6b) 

Z X C 

a = — = (— sin 0, cos 0, 0) , (6c) 

sin 9 

b = c x a = (— cos 9 cos 0, — cos 9 sin 0, sin 9) . (6d) 

As a notation to be used throughout this work, subscripts z, c, a or b on a vector indicate 
its orthogonal projection on the z, c, a or b directions and subscript c indicates the vector's 
projection on the plane perpendicular to c . 

The particle's angular velocity vector u) may be decomposed into a sum of two vectors, 
perpendicular and parallel to c, respectively, 

uj = U)c + uj c c , 

with 

LJc = c x c = (—9 sin — sin 9 cos 9 cos 0, 9 cos — sin 9 cos 9 sin 0, sin 2 0) 

and 

Cg> c = "0 + COS 9 . 

The orthogonal projection of <jJ 5 on the z-axis is 

ujcz = w s • z = sin 2 # , 

and the orthogonal projection of a; (or of u E ) on the direction perpendicular to the cz plane 
is 

u a = u) • a = <jj 5 - o, = 9 . 

Thus we see that the dissipative torques present in Eqs. fl5a|), ( |5l5| ) and fl5c|) are given by u a , 
ujcz and uj c , respectively, times the dissipation parameters A or A'. 

The noise torques will be treated along these same lines. We start by defining the noise 
torque vector by its orthogonal components, 

T = T a a + T b b + T c c . 

The noise becomes completely defined by stating the statistics of its three components. 
The usual procedure is to consider them as statistically independent, Gaussian white noise. 
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This is, however, not a necessary assumption and we leave it open for future modeling. 
What we need now is to know how the three components come into Eqs. (H). Guided by 
the above decomposition of the dissipative torque, we are led to identify 

r<j = r a , 

rv = r gz = r g • z = r b sin o , 
rv = r c . 

Before we proceed to deduce the equations of motion for the general case of magnetic 
particles in ferrofluids we show, in the next section, how to obtain, from Eqs. (||), the 
equations of motion for the spherical coordinates of a mono-domain magnetic moment. 



III. EQUATIONS OF MOTION FOR A MAGNETIC 

MOMENT 

The magnetic moment /x of a mono-domain particle is related to its internal angular 
momentum S by fi = 7S, where 7 is the gyromagnetic factor. Although the modulus 
S of S is taken as constant in most works on superparamagnetism and magnetic fluids, 
for very small particles its oscillation may be significant and we prefer to allow it to be 
time dependent. The modern technology allows the preparation of samples with magnetic 
particles whose diameters are smaller than 20A |H| and superparamagnetic clusters con- 



taining only 12 magnetic atoms have also been reported [TJ|. We can model the magnetic 
moment by a rotating charged particle, in the limit of zero moments of inertia, l\ — ► 0, 
I3 — ► 0, and ip — > 00 so that I 3 ip = S. Because in the next section we will work with 
the joint system, a particle and its fluctuating magnetic moment, we write the general- 
ized coordinates, potential energy, dissipative and noise torques, with a notation distinct 
from that corresponding to the particle. Namely, we make the following substitutions: 
6 -> <f> -»• <p, I 3 ip -> S, V -»• W, A ^ £, A' -> f ' and T -> T. We also introduce two 
modifications in the equation corresponding to Eq. fl5c|) , namely, we write S — S instead of 
S in the dissipative term and introduce a torque W s , whose origin will be explained below. 
In the said limit and with the new notation the system of Eqs. @ becomes: 

S (p sintf + £Q + We = % , (7a) 
Scosti- S $sm$ + £^sm 2 $ + W v = % , (7b) 
S + aS-S ) + W s = T s . (7c) 



Here we have written S — So, instead of S, in the dissipation term of Eq. ( |7cD to account for 
the fact that the relaxation of the fluctuations of S is towards a most probable (equilibrium) 
value So and not towards 0. It may appear strange that, even though we have derived the 
equations of motion for S from the equations of motion for a symmetric particle, in a con- 
venient limit, we have now to add a term "ad hoc" (So), which does not have an equivalent 
in the particle's equations. This is so because in classical physics the equilibrium magne- 
tization is always zero. Non-zero equilibrium magnetic moments can only exist because of 
the quantum mechanical nature of matter and, therefore, cannot be deduced from a pure 
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classical approach. The torque W s was introduced because a crystal field may have an ef- 
fective interaction with /j,, with origin in an orbital contribution to S [24], with a possible 
torque component parallel to S. There is not an equivalent term in Eq. ( pcD because of the 
assumed axial symmetry of the particle. 

It is interesting to study the behavior of Eqs. ([/]) in the absence of noise, % = and with 
W s = 0. Eq. ([TcD has then the trivial stationary solution S = Sq. Assuming this constant 
value for 5* in Eqs. ( |7a| ) and flTE| ) they reduce to 



SoVsm<d + £<d + W# = Q , 

-S $sm$ + £ip sin 2 ■& + W v = . 



(8a) 
(8b) 



The conservative torques, — W$ and —W v , have, usually, contributions from two different 
origins, the interaction of S with a crystalline, anisotropy field and/or with a magnetic field, 
which can also be of several different origins. In the case of magnetic field, H, the potential 
energy is W = — n ■ H. With a little of algebraic work one can show, in this case, that the 
set of Eqs. (H) is equivalent to the well known Gilbert's equation M, 



~dt 



7 fi x 



H 



£ d/j, 

/i 2 dt 



for n = and S = So- This equation was used by W. F. Brown JTT[] as a starting point 
for his stochastic theory of superparamagnetism, where he assumed the magnetic field H 
to contain a noise term. A more general theory for superparamagnetism, which allows also 
for oscillations on the modulus /i = 'yS of the magnetic moment, was worked out by Ricci 
and Scherer fT7HT9!, based on the set of Eqs. (^). For this reason we will not continue to 



explore the consequences of Eqs. (0) in the present paper, turning, instead, to the more 
general ferrofluid, where the rotations of the mechanical particle are taken into account, in 
addition to the motion of S relative to the particle. 



IV. THE GENERAL FERROFLUID 



In recent years several researchers JT]fJ,|T^,^] have drawn attention to the importance 



of the motion of the magnetic particle, its inertia and viscous interaction with the fluid, to 
the dynamic magnetic susceptibility of ferrofluids. A theoretical treatment of this problem, 
which is both, more fundamental and more general than those previously published, follows 
naturally from the context described above. 

Taken together, the systems of Eqs. @ and (0) contain all the degrees of freedom 
relevant to the problem. To the potential energy terms, V in Eqs. ([5]) and W in Eqs. (0), 
the interaction energy between the magnetic moment and the particle, which we will denote 
by U, has to be added. Due to the particle's symmetry, this term can only depend on S 
and on the angle between S and the symmetry axis, c. It is convenient to define another 
orthogonal set of unit vectors, related to the direction of the magnetic moment, namely, s, 
in the S direction, u, perpendicular to the sz-plane and v, perpendicular to the su-plane, 
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s 

s = — = (sin $ cos 99, sin 1? sin 99, cos #) , (9a) 

2 X S . . / -i \ 

u = — = (— sin 09, cos 09, 0) , (9b) 

sin 17 

v = s x u = (— cos i9 cos 09, — cos sin ip, sin -(9) . (9c) 

The interaction energy U can then be written as U(S, s ■ c). In principle the particle can 
interact also with other fields, besides H, as is the case if it has an electric dipole and an 
electric field is present. For this reason we keep also the potential energy V(6, <p) in the new 
set of equations. 

The dissipative interaction associated with the rotation of S relative to the particle will 
be written in terms of the relative angular velocity vector. Since only rotations perpendicular 
to S can lead to a meaningful interaction torque with origin on the relative motion, we define 
the relative angular velocity uj t as 

where 

VJ = S X s 

is the angular velocity of rotation of the magnetic moment with respect to the laboratory 
and 

iVs = SXLJXS = <jL! — (s ■ uj)s . 

is the orthogonal projection of the particle's angular velocity uj on the plane perpendicular to 
S. The dissipative interaction torque on the particle is then +£ oj r . The plus sign is because 
of the way we defined cj r , where the particle's angular velocity appears with a minus sign. 
Guided by the interpretation of the dissipative torque terms of Eqs. @ in terms of angular 
velocity components, as explained bellow the said equations, we write down immediately 
the dissipative torque terms to be added to the left-hand sides (therefore, with a — sign) of 
Eqs. (^), namely 

-£ u ra = -£ u) r ■ a , 

-f ^rcz = -f K - (UV • C) C] • Z = -f {U) rz - U rc COS 6 , 

-£ uj tc = — £ uj t ■ c . 

Of course, all this scalar products, as well as those which follow, in the next equations, may 
be easily written as functions of the four angles 8, <fr, •& and 09 and their time derivatives, 
by using Eqs. @ and @. However, because scalar products are very easily handled in 
numerical procedures, we prefer to leave them in this form. 

Clearly, the torque on the magnetic moment, due to the relative motion, is the "reaction" 
to the torque on the particle, i.e., it is equal to — £u> n and, in place of £ $ and £ 09 sin 2 1? in 
Eqs. ([?D we shall use (remembering that 

£ OJru = £ W r ■ U , 

£ u rz = £ u; r ■ z . 
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No term coming from the relative angular velocity u) r has to be added to Eq. (|7c ) because 
uj r is perpendicular to S. However, there is the term £' (S — Sq) already present in that 
equation, with origin in the (quantum) fluctuations of S, and this term will be kept. Since 
angular momentum has to be conserved, its reaction counterpart on the particle has to be 
added to Eqs. (|). Calling 

n = (s - So) s , 

the terms to be added to the left-hand sides of Eqs. (|5p are 
1l a = -en-a = (S-S )s-a, 

n- cz = -e [k-(k-c)c]-z = -r (s - s ) [ s - ( S ■ C ) C ] • % , 
-e' n c =-en-c= -e (s-s )s-c. 

The noise torques of interaction between the particle and the magnetic moment can be 
written down along the same lines of procedure as done for the noise torques of the fluid 
on the particle, at the end of section II. We assume three orthogonal, independent, noise 
torque vectors, along the unit vectors defined with respect to the direction of the magnetic 
moment: 

T = T S s + T u u + %v . (10) 

Being T the torque on the magnetic moment, then the torque on the particle is —T. 
Following the same line of reasoning as done before, we identify the torques in Eqs. (0): 

% = % z = % sin i? , 

z = z . 

Correspondingly, we the following terms have to be added to the right-hand-sides of Eqs. 
(I): 

Z = -Z = -T-a, 

Zp = -T Ez = —[T — (T ■ c) c) ■ z , 

Zp = -Z = -Z ■ c . 

Therefore, the state of the composed system, the particle and its magnetic moment, is 
described by the 6 generalized coordinates, 9, <fi, '0,1?, y 9 and S, whose dynamical behavior 
is governed by the following set of coupled differential equations: 
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h(0 - 4> 2 sin6>cos6») + I 3 4> (ip + <ficos9) sin9 + A - & ra + 

-en a + Ve + U e = T a -T a , (11a) 
I 1 ((f)sm 2 9 + 2 <p 6 sin 6 cos 6) + I 3 cos0 ^-(ip + <j)cos9) + 

-I 3 (ip + cos 0)0 sin + A sin 2 - fu; ra - f ' 7L- CZ + V < f, + U 4> = V b sin 0-%, , (lib) 

^(V* + 0cos#) + A' (V> + 0cos#) - £u; rc - £' ^ c = T c - T c , (11c) 
5 ip sin + £w rM + W# + Uo = +T U , (lid) 
S'cos^ - S i9sini? + ^ rgz + W 7 ^ + t/" v = +T 5z , (lie) 

S + g{S-S ) + U s = T a . (llf) 

This set of six equations is of very general applicability on ferrofluids. It allows for a large 
variety of modeling: There are three independent conservative interaction potentials, V, U 
and W, four dissipative parameters, A, A', £, and and also the noise torques T and 
7", whose statistical properties are open for modeling. Particle-particle interaction was not 
explicitly taken into account, but on a mean-field approximation it can be included in V 
and/or in W. 

As we mentioned before, in most cases of practical interest the fluctuations in the modulus 
S can be neglected. In this case Eqs. ( |TTD become simpler, in several respects: Eq. ( |llt| ) 
ceases to exist, all terms in S and in £' become zero and the noise term % in Eq. (|T0|) and 
its contribution in Eqs. ([TT| ) also vanish. 

Two interesting limit situations are readily obtained from Eqs(|ll|), the "superparamag- 
netic" limit, for which the particle's coordinates, 0, <f), and ip, are taken as constants, so 
that the system reduces to the last three equations, and the "blocked" limit (also called 
"Brownian" limit 0] or "inertial limit" [§), when the magnetic moment is blocked along 



the particle's symmetry direction, i.e., d = 9 and ip = 0, but the particle can rotate inside 
the fluid. The superparamagnetic limit has been treated in three previous papers by Ricci 
and Scherer JT^-IT^] and also by other authors. In the remaining sections of this paper we 



will deal with the "blocked" limit. 



V. DYNAMICS OF THE MAGNETIC MOMENT IN THE BLOCKED LIMIT 

We consider now the situation in which the magnetic moment is blocked along the par- 
ticle's symmetry axis. This may happen because the sample is kept below the "blocking 
temperature" Tb p2| or because the material is so highly anisotropic that the magnetic 



moments only exists parallel to the easy axis [[[(J. The particle is still immersed in a fluid 



carrier, being able to rotate, together with its magnetic moment. 
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In terms of the set of Eqs. ([0]), the blocked limit is obtained by assuming an interaction 
potential U of the form — Uq5(s — c), with Uq — > oo, so that the only states energetically 
possible are those with s = c, i.e., d = 9 and ip = 0. By summing Eq. ( |lla| ) with Eq. ( |Hd| ) 
and Eq. (|llb|) with Eq. ( p.le| ) the interaction terms Uq and U$ as well as and U v cancel 
out. The terms containing u ra , w r5z , 7£ a , T^-cz, T a , and T gz become identically zero, and 
TZ C becomes (S — So)- Choosing and <fi to denote the common polar angles, the system of 
equations becomes: 

h(9 - 4> 2 sin 9 cos 9) + J 3 <fi (?/> + 0cos0) sin0 + A 9 + V e + 

+S0sin0 + W e = T a , (12a) 

I 1 ((f)sm 2 9 + 2 <j)9sm9cos9) + I 3 cos9 ^-(ip + 0cos6») + 

at 

-I 3 (tp + cos 0)5 sin + A 0sin 2 + V^ + ScosO + 

-5 0sin0 + ^ = r b sin0, (12b) 

hj t $ + 0cos0) + A' (V> + 0cos0) - i\S - S ) = T C -T C} (12c) 

S + e(S-S ) = +T c . (12d) 

We will not explore, in the present paper, the system of Eqs. (0) in all its generality. In 
the following we will consider only the cases when the noise terms do not need to be taken 
into account explicitly. The explicit presence of white noise torques makes of Eqs. (|TT|) and 



(12]) Ito-Langevin systems and their treatment demands the use of stochastic calculus. This 
will be treated in a future paper. Several circumstances may be devised where neglecting 
the noise is a good approximation. One of them is when the system is drawn far from 
equilibrium, in presence of a strong magnetic field, with fiH ^> fe^T. Then the relaxation 
to the new equilibrium state, with /x approximately parallel to H, proceeds without an 
important influence of the noise. The Bloch's equations of magnetic resonance are based 
on this idea: they contain relaxation terms (with relaxation times T\ and T 2 ), but no noise 
terms. Another interesting circumstance is in the context of linear response theory. The 
usual formulation considers the noise sufficiently weak for its effect to be only in establishing 
an initial thermal equilibrium. The perturbing field is then introduced adiabatically, i.e., 
with the system disconnected from the thermal bath. A similar approach to linear response, 
however based on the equations of motion, instead of based on the Liouville equation for 
the probability distribution, which is the case of Kubo's linear response theory PBj, will be 



presented in the next section. 

In the following we will also assume a constant modulus for the magnetic moment, i.e., 
S = So, and for the interaction potential we consider only W = — /Lt ■ H = —^SoS ■ H. For 
simplicity, only a constant field, H = H z, will be considered now, but in the section on 
linear response we will introduce also a time dependent transversal field. 

With this simplifications, the system of Eqs. (|12|) becomes 
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h{d - <p 2 sin0cos0) + h (j> (ip + 0cos0) sin0 + A + 
+5*0 sin 9 + ■jSqHo sin = , 



(13a) 



/i(0sin 2 + 2 sin cos 0) + I 3 cos9 ^-(ip + 0cos 
-I 3 (ip + <j) cos 0)0 sin + A sin 2 - S 9 sin = 



h-rid* + 0cos0) + A' (V> + 0cos0) = 



(13b) 
(13c) 



Eq. (13c) shows that, under the circumstances considered, and for any value that the 
function if) + 0cos0 may have, due to initial conditions, it will relax exponentially to zero. 
Since we will consider, in what follows, only stationary initial conditions, its value will be 
taken as identically zero. This simplifies the above set of equations to 



h( 



sin cos 0) + So 4> sin + A 6 1 = —^SqHq sin 



J 1 (0sin 2 + 2 0sin0cos0) - S sin + A 0sin 2 = . 



(14a) 
(14b) 



This equations may be solved numerically, for arbitrarily given parameters, h, Sq, A and 
jH and arbitrary initial conditions, by using standard methods. An example of solution 
is shown in Fig.l, in form of a rotational trajectory of the magnetic moment, drawn over a 
sphere to help visualization. The origin of the magnetic moment vector is at the center of the 
sphere and its head follows the trajectory on the surface, the magnetic field H is parallel 
to the line from the south to the north pole. For Fig.l-a the initial velocities vq = 9(to) and 
wo = 4>{to) were arbitrarily chosen, for Figs.l-b and 1-c the initial velocities were calculated 
from Eqs. ([37]) and (|38|) of Sec. VIII. The dissipation parameter A for Fig. 1-c is 100 times 
the value used for Figs.l-a and 1-b. All other parameters and also the total time interval 
are the same for the three figures. 




Fig.l: Trajectory of the head of the magnetic moment vector in its rotation around the constant 
magnetic field (see text). 1-a: for arbitrarily chosen initial velocities; 1-b: for initial velocities given 
according to Sec. VIII; 1-c: same as 1-b, but with dissipation parameter, A, 100 times larger. 
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VI. A SIMPLE APPROACH TO LINEAR RESPONSE 

THEORY 



The standard Linear Response Theory p3| is based on the classical or quantum Liouville 
Equation, for the probability distribution function or the density matrix, respectively, for 
systems in absence of noise, or based on the Fokker-Planck Equation in the case of stochastic 
processes [19]. A simpler approach, which has some limitations, but is appropriate for the 
purposes of this work, based on the equations of motion of the system, is as follows. 



A. The Linear Equations for the Perturbation 

Let us assume that the system's coordinates satisfy the equation 

Q(t) = P(Q,t), (15) 

where Q = (gi, q2, ■ • ■ , q n ) and where P = (Pi, P2, • ■ ■ , P n ) are functions of the coordinates 
and of the time. We assume further that the only explicit time dependence of P comes from 
an applied perturbing field, 

F(t) = (Fi(t), • ■ ■ , F m (t)) , 

in the form 

P(Q,t) = P°(Q)+A(Q).F(t), (16) 

where A is an n x m matrix. The perturbation F(t) is assumed sufficiently weak for its 
effect in deviating Q(t) from its unperturbed values to be well approximated by a a linear 
functional. This approach is exact in the limit F — > 0, which defines the initial susceptibility. 
In this equation and in what follows an upper index (like in P ) will indicate "unperturbed 
values", while a lower index (like in Q ) will be used for initial values. Therefore Q (t) 
indicates the solution of the unperturbed equation, with given initial conditions, 

Q°(t) = P°(Q°) , Q°(t ) = Q° = (?oi, • • • , Qon) ■ (17) 



The solution of Eq. (|T5|) will be written as 

Q(t) = Q°(t) + X(t), (18) 



where X(t) is a linear functional of F. Substituting Eq. (18) into Eq. (|15D , using Eq. 
and keeping only linear terms in F follows 

X(t) = K{Q\t)) ■ X(t) + A(Q°(t)) ■ F(t) , (19) 

where the matrix elements of K are 

dP° 

*» - 1± ■ < 20 > 
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Eq. ( |I~9"D is a linear, inhomogeneous, system of equations for X, with time dependent 
coefficients K^t) = K^j(Q (t)). The corresponding homogeneous equation, 

Y = K-Y (21) 

has the solution 

Y(t)=M(t,t )-Y , 
where the matrix M is formally given by 

M(t, t ) = exp( f K(t') dt') . (22) 



The general solution of Eq. ([19]) may be written formally as 

X(t) = M(t, t ) ■ X + I M(t, t') ■ A(t') ■ F(t') dt 1 . 

Jt 

Since X(t) has to be a linear functional of F, it follows that X = 0. Thus 

X(t) = f M(t, t') ■ A(t') ■ F(t') dt' . (23) 



In numerical procedures, it is often simpler to solve Eq. (^T|) then to work with Eq. 
to obtain the matrix elements of M. To obtain Mij(t,to) from the solutions of Eq. (^) we 
define the set of n unit vectors of dimensionality n, 

rS = (i,o,...,o) 

Yl = (0, 1,-", 0) 



n = (o,o,---,i) 

The solution of Eq. ([21]), with initial condition Y(to) = Y will be denoted by Y l (t), 



i.e., 



Y i (t) = M(t,t )-Y . 

Therefore, the j th component of Y l (t) is 

Y>{t) = Mfifato) , (24) 

from what follows that to obtain the matrix elements Mji(t,to) one solves the linear set of 
Eqs. ([21]) with the Y as initial vectors. 
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B. Response Function and Susceptibility Matrices 

Consider the dynamical variable (observable) B(Q) = (Bi, ■ • ■ , B m ). Its ensemble aver- 
age over the initial equilibrium distribution will be denoted by (B(Q)) Q . The "response" of 
B to the perturbing field F(t) is defined by 

5B(t) = (B(Q(t)))o ~ (B(Q°(t))) Q . 
Expanding the first term above around Q° to first order in X and using Eq. (|23|) follows 

SB(t) = (vB(t) ■ f M(t, f) ■ A(t') ■ F{t') dt'\ , (25) 

\ J to I 

where V-B is an m x n matrix, with elements 

— -) , fc=l..-m, i = l---n. (26) 

°1i J Q=Q0(t) 

The "response function matrix" (m x m) $ is defined by 

5B(t) = f $(t-t) ■ F(t')dt' . (27) 

Jto 

Therefore, comparing Eqs. (]27f) and (^), we identify 

^t-t') = (VB(t)-M(t,t').A(t')) o , (28) 

which is function of t — t' and not of t and t', independently, because the average is over 
the equilibrium distribution, for which absolute times are meaningless. Therefore, we can 
choose t' = to = in Eq. (p8|) and write the matrix elements of $(£) as 

= J] (B^Mijit, 0)A jl (0)) . (29) 
The complex susceptibility is defined as the Fourier-Laplace transform of $, 

("OO 

Xw(^) = Km / exp(i^t / - et')$ k i(t') dt' . (30) 

e^0+ JO 

In the next section we apply the concepts and results presented above to the ferrofluid 
in the blocked limit. 



VII. DYNAMICAL SUSCEPTIBILITY OF THE BLOCKED 

FERROFLUID 

The starting point to apply the linear response approach of last section is the set of 
equations of motion for the system. We rewrite Eqs. fll4]) in the form of Eq. (|I5]), by 
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defining the variables v = 9 and w = <ft. To simplify the notation we also introduce S 
So/ii, X — X/Ix and write Eqs. ( |Tjj ) in the form 



e -. 

0: 



w 2 sin 9 cos 6* — S w sin — A v — 5 7 H sin 
= —2 v w cot + S v esc 6* — A w . 



(31a) 
(31b) 
(31c) 
(31d) 



By comparison of Eqs. ([H]) with Eq. ([T7|) we see that the components of the vector 



P (Q ) are just the RHS's of Eqs. fl3~T|). From Eq. (|20|) we obtain the matrix elements 
Kij(t), i,j = 9, 0, v, w. The only non zero elements are 

K dv = 1 , 

K v e 
K 



K 

K 
K 



w (1 — 2 sin 
-A, 

: (2 w cos 9 — S) sin 9 , 
■ (2vw — S v cos 6 1 ) / sin 2 9 

- (-2wcosfl + S)/sinfl , 

- —2 1> cot 9 — A , 



5(7 i^o + cos 6* 



where 9 = 9(t), v = v(t) and w = w(t) are the solutions of Eqs. For any given set of 

initial values 9o, vq and wo it corresponds a set of functions Kij(t) (independent of 0o) and, 
from Eqs, fl2"ID and (|4]) follows the corresponding Y/ (t) and My(£). 

We assume now that a time-dependent perturbing magnetic field is applied perpendicular 
to the z-axis, i.e., H(t) = (H x (t), H y (t), 0). The interaction energy of this field with the 
particles magnetic moment is 

W = —(/, ■ H = —S 7 H x sin 9 cos — So 7 H y sin 9 sin . 

The terms to be added to Eqs. (|31cj ) and (|31d|) are 



-ldW 
~h~d9 
-1 dW 



•yS cos9(H x cos <fi + H y sin < 
7 S 



Ji sin^ 



sin 9 



-H x sin + H y cos 1 



By adding this terms to the RHS of Eqs. ( |31c|) and (|31d|) , respectively, comparing with Eqs. 
(|i5|) and (16), and identifying F with H, we can write down the matrix A: 



A 



( A-Q X A0y ^ 

A i 4 j 

■^(px -^fpy 

A A 

V ^tra A wy / 



( 






7S cos 6* cos 7S cos 6* sin 

V — 7>S'sin 0/sin# yS cos(f)/ sin9 J 



(32) 



15 



To calculate the response functions <&ki{t) by Eq. ( p9|) only the matrix elements at time zero 
are needed. For this purpose we substitute 9 and in Eq. ( 32|) by 9 and O . 

To complete the argument of the average in Eq. ( p9| ) we still need the B^^t). As the 
observable vector B we choose the projection of the magnetic moment /x on the xy-plane, 
i.e., B = (B x , By), with 

B x = 7 So sin 6* cos , 
By = j S sin # sin . 



From Eq. (p6l) we then get 



B x ,e = 7 S cos # cos , 
B x ^ = -7 So sin 9 sin , 
-£>j/,e = 7 S cos 9 sin , 
-By,^ = 7 So sin # cos , 

and all other B^i are zero. 

For the response function <f> xx (t) we obtain, from Eq. 

®xx(t) — (B Xt gMg v A vx + B x fiM ew A wx + B X ^M^ V A VX + B x ^M ( j )W A wx ) 



(33a) 
(33b) 
(33c) 
(33d) 



o ' 



(34) 



with time arguments B Xji (t), Mij(t,0) and Aj X (0). Substituting the B xi and Aj X by their 
values as given by Eqs. (|33| ) and fl32|), writing <fi(t) = 0o + A0, where A0 is independent 
of 0o (because the unperturbed equations do not contain 0) and using some trigonometric 
relations, the average over 0o in Eq. (B3) may be done analytically, resulting in 



'j 2 S S 



cos 9 

cos 9 cos 6*o cos A0 Mq v H — : — — si 1 1 



sin 9q 



M ew + 



sin 9 cos 9 sin A0 Ma v + 



sin# 
sin 9n 



cos A0 M t 



dill! 



(35) 



J / o 



By the same procedure we obtain the other response functions: 



7 2 S S 



cos 9 

cos 9 cos 6*o cos A0 M dv H — : — — sin A0 M dw + 



— sin 9 cos 9q sin A0 M^ v + 



sin 6*o 
sin 9 



sin 9q 



cos A0 M« 



(36) 



Eq. 



also shows that Onsager's relation 

% x (t,-H ) = <S> xy (t,H ) , 

—H . We also obtain <3> TO 



is satisfied, because $ xy is odd under H 
obvious result, from symmetry. 



$xx, which is an 
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VIII. NUMERICAL RESULTS AND CONCLUSIONS 



The equilibrium average indicated in Eqs. (pq ) and fl3"E| ) are, in principle, to be done 
over 9 , v , w and all particle's independent parameters (polydispersity). The average over 
0o was already performed analytically. However, the deviations of v and w from their most 
probable values are rapid fluctuations due to the Brownian noise, which is not present in 
our system of equations. Therefore, the choice of initial distribution, to be consistent with 
our approximation of neglecting the noise, should not include this fluctuations. We will use 
Boltzmann's equilibrium distribution for 6*0 and, for every selected value of 9 , we chose v 
and w from an approximate solution of Eqs. ( |31cp and ( ^ld| ) around t = 0, calculated in 
the following way: 

1) Assume that if v is not zero at t — 0, when we disconnect the noise, it relaxes to zero 
according to the equation v = —A v. Eq. Q31c| ) then leads to 



S-^S 2 + 4 7 tf ,Scos£ 

= = 7, , 37 

2 cos # 

where we have chosen the — sign because wq has to vanish for H = 0. 

2) Analogously, we assume that if w is not Wq, it relaxes to Wq according to the equation 
w = —A (w — w ). Eq. ( pi d| ) then leads to 

v = - . (38) 

S — 2 w cosOq 

This prescription was used in Sec.V for the initial velocities in the calculation which led to 
Figs. 1-b and 1-c. We remark that this choice of Wq is meaningful only if 

5 + 4 7 #ocosfl > . (39) 

However, if 7 So Ho/hsT is not too small, the Boltzmann distribution becomes negligible for 
values of 9 Q such that Eq. Q3DD is not satisfied. 

It is also important to examine under which circumstances the approximation made in Eq. 
(|l3|), and in all that follows that equation, the neglecting of the noise terms, is appropriate. 
We are specially interested in obtaining the dynamical susceptibility of the system, and 
therefore we will examine the consequences of that approximation in the context of linear 
response. Two characteristic times are of importance: The longitudinal relaxation time 
T\ pa where the average is over all magnetic particles, and the transverse relaxation 

time T 2 , which is the time taken by the response functions to become approximately zero. 
We have borrowed the notation 7\ and T 2 from Magnetic Resonance (MR) because of the 
similarity of their meanings here and in NMR or EPR. The longitudinal (parallel to Hq) 
relaxation, called "spin- lattice relaxation" in MR, occurs in a time Ti, via energy loss to 
the environment, due to the dissipative torque. The transverse relaxation, characterized by 
the vanishing of <& xx (t) and called "spin-spin relaxation" in MR, occurs in a time T 2 , due to 
the loss of phase coherence between the responses of the different particles to a pulse of the 
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perturbing field at t = 0, in their rotations around H . Since we assume that the system 
is very close to thermal equilibrium, at a given temperature, and since the longitudinal 
relaxation takes the system out of the initial equilibrium distribution corresponding to that 
temperature, just because the thermal noise is neglected in the equations of motion, our 
calculation of the response functions and susceptibility is only a good approximation if 
T 2 < T x . 

There are several sources of transverse relaxation. One of them is, of course, the noise, 
which is being neglected in this approximation. The different initial angles 8 and different 
particle's parameters Ji, A, and So (polydispersity) lead to different frequencies (see Eqs. 
(|3T|) and (|38|)) and, consequently, to loss of phase coherence and to transverse relaxation. 
These are taken into account in the averaging procedure on Eqs. ([35]) and (|36|) . 

To obtain the functions 0(t), A<p(t) and Mij(t) we have to solve the systems of differential 
equations, Eqs. (^Tj) and (|3l"l) . We have done it with the Runge-Kutta algorithm, in a work- 
station. The particle's parameters, field intensity and time unit have been arbitrarily chosen 
and kept the same in all calculations whose results are shown in the figures, except where 
explicitly stated. 

Polydispersity was treated for particles made of the same material and having the same 
shape, assuming a uniform distribution of a linear dimension, r, in an interval of size Ar, 
i.e., r was chosen to be uniformly distributed in the interval (1 — |Ar, 1 + |Ar). The other 
particle's parameters were scaled accordingly, namely, 

So oc r 3 , l\ oc r 5 , 

S = Sq/Ii oc r~ 2 , 

A oc r 3 , A = A/ii oc r~ 2 . 

The average over the initial angle 8 Q was weighted with the Boltzmann equilibrium 
distribution, 



P(6> ) oc sin^o exp(So7-ffo cos^o/Zc^T) , 

5, so that P( 



and the temperature was chosen so that So 7 H cos 6 /k B T 
mum around 71/6. 



has a maxi- 




Fig.2: & xx [t) versus time. Line a: Ar = 
Line b: Ar = 0.05; Line c: Ar = 0.20. 



Fig. 2 shows the diagonal response func- 
tion <fr xx (t) versus time, for three differ- 
ent polydispersity ranges, Ar = 0, 0.05 
and 0.20, for the curves a, b and 
c, respectively. This figure confirms 
what we said above that polydisper- 
sity causes the relaxation time T 2 to 
decrease. Their values, for the curves 
a, b and c, may be estimated to be ap- 
proximately T 2 ~ 0.2, 0.1 and 0.05, re- 
spectively. The longitudinal relaxation 
time Ti cannot be estimated from this 
curves, but has to be calculated to- 
gether with the numerical solution of 
Eqs. (|3Tp, by using 7\ w (0)- 1 . 
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The result is not strongly dependent on polydispersity, giving, in the present case, T\ m 
20. Therefore the condition stated above for the appropriateness of the approximation of 
neglecting the explicit presence of noise in the equations of motion, T 2 <C T l5 is amply 
satisfied for the parameters used here. 

Fig. 3 and Fig.4 show the real and imaginary parts of the susceptibility respectively, 
corresponding to the response functions of Fig. 2. Among other information, Fig.4 shows 
clearly that the broader the dispersity of particle's size, the broader also the resonance line, 
as one should expect. 




200 400 

frequency (0 



600 



Fig. 3: Real part of the susceptybility 
x(w), versus u, corresponding to the re- 
sponse functions of Fig.2. 



Si 




200 400 

frequency co 



600 



Fig.4: Imaginary part of the susceptybil- 
ity x( w )> versus u, corresponding to the 
response functions of Fig.2. 




0.0 .05 .10 .15 .20 

time 

Fig.5: Solid line: <& xy (t); dot-dash line: 
®xx(t), the same as in Fig.2-a. 




300 400 



frequency co 

Fig.6: Imaginary part of the suceptibility 
for different values of the moment of in- 
ercia, I\. Full line: l\ = 1.0; dashed line: 
h = 0.25; dot-dash line: h = 0.10; dotted 
line: h = 0.05. 
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Fig. 5 shows the non diagonal response function $ X y(t) (solid line), for a monodisperse 
sample, plotted together with <& xx (t) (dot-dash line) for comparison. In Fig. 6 we compare the 
resonance frequency (center of the resonance peak in the imaginary part of the susceptibility), 
for different values of the moment of inertia, I±, keeping constant all other parameters. 
The lowest value, I\ = 0.05 (dotted line) is the same as used in the previous figures, for 
monodisperse samples. The other curves correspond to l\ = 0.10 (dot-dash line), 0.25 
(dashed line) and 1.0 (full line). Clearly, the heavier the particles, the lower the resonance 
frequency. 
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